import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import minimize
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
import scipy
from sqlalchemy import create_engine, inspect
from sqlalchemy.orm import sessionmaker
import datetime
Para la base de datos se utiliza el servicio de https://neon.tech/, el cual se utiliza por integrar PostGIS sin costo y ofrecer una capacidad suficiente para pruebas.
db_url = "postgresql://test_owner:lpPrcHty8i9V@ep-royal-glitter-a50h4o6a.us-east-2.aws.neon.tech/test?sslmode=require"
engine = create_engine(db_url)
Se guardó toda la información de la Especie con identificador 1 para hacer pruebas.
query = "SELECT * FROM data_especie1;"
df = pd.read_sql_query(query, engine)
df
| PARCELA # | ROTACIÓN # | ZONA | AÑO PLANTACIÓN | EDAD | MEDICIÓN # | FECHA MEDICIÓN | ESPECIE | PROCEDENCIA GENÉTICA | ARBOL # | ... | DAP MEDIO [cm] | AREA BASAL [m2/ha] | H TOTAL [m] | H DOMINANTE [m] | H PODA [m] | N/ha PODADOS [#] | VOLUMEN EN PIÉ d.u. = 5 cm. [m3/ha scc]. | IMA s/raleo [m3/ha/año]. | Observación | Clasificacion Densidad | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1004 | 2012 | 6.043836 | 2 | 2018-06-20 00:00:00 | 1 | Livingston HSCA T416 | 1 | ... | 17.229788 | 14.112918 | 9.768085 | 10.925 | 3.555319 | None | 51.132477 | 8.460269 | None | 3 |
| 1 | 1 | 1 | 1004 | 2012 | 6.043836 | 2 | 2018-06-20 00:00:00 | 1 | Livingston HSCA T416 | 2 | ... | 17.229788 | 14.112918 | 9.768085 | 10.925 | 3.555319 | None | 51.132477 | 8.460269 | None | 3 |
| 2 | 1 | 1 | 1004 | 2012 | 6.043836 | 2 | 2018-06-20 00:00:00 | 1 | Livingston HSCA T416 | 3 | ... | 17.229788 | 14.112918 | 9.768085 | 10.925 | 3.555319 | None | 51.132477 | 8.460269 | None | 3 |
| 3 | 1 | 1 | 1004 | 2012 | 6.043836 | 2 | 2018-06-20 00:00:00 | 1 | Livingston HSCA T416 | 4 | ... | 17.229788 | 14.112918 | 9.768085 | 10.925 | 3.555319 | None | 51.132477 | 8.460269 | None | 3 |
| 4 | 1 | 1 | 1004 | 2012 | 6.043836 | 2 | 2018-06-20 00:00:00 | 1 | Livingston HSCA T416 | 5 | ... | 17.229788 | 14.112918 | 9.768085 | 10.925 | 3.555319 | None | 51.132477 | 8.460269 | None | 3 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 383240 | 106 | 1 | 1004 | 2012 | 10.909589 | 7 | 2023-07-10 00:00:00 | 1 | HSCB R | 35 | ... | 28.300000 | 25.390000 | 17.600000 | 18.000 | 4.900000 | None | 173.100000 | 15.866775 | None | 2 |
| 383241 | 106 | 1 | 1004 | 2012 | 10.909589 | 7 | 2023-07-10 00:00:00 | 1 | HSCB R | 36 | ... | 28.300000 | 25.390000 | 17.600000 | 18.000 | 4.900000 | None | 173.100000 | 15.866775 | None | 2 |
| 383242 | 106 | 1 | 1004 | 2012 | 10.909589 | 7 | 2023-07-10 00:00:00 | 1 | HSCB R | 37 | ... | 28.300000 | 25.390000 | 17.600000 | 18.000 | 4.900000 | None | 173.100000 | 15.866775 | None | 2 |
| 383243 | 106 | 1 | 1004 | 2012 | 10.909589 | 7 | 2023-07-10 00:00:00 | 1 | HSCB R | 38 | ... | 28.300000 | 25.390000 | 17.600000 | 18.000 | 4.900000 | None | 173.100000 | 15.866775 | None | 2 |
| 383244 | 106 | 1 | 1004 | 2012 | 10.909589 | 7 | 2023-07-10 00:00:00 | 1 | HSCB R | 39 | ... | 28.300000 | 25.390000 | 17.600000 | 18.000 | 4.900000 | None | 173.100000 | 15.866775 | None | 2 |
383245 rows × 42 columns
Para el crecimiento se definieron 12 modelos de crecimiento. El ejercicio consiste en encontrar, a partir de los datos, los mejores valores constantes para cada modelo.
#Sweda
#N(t) = N0 * e^(r(t-t0))
def SWEDA(x,N0,r,t0):
return N0 * np.exp(r*(x - t0))
#return x * N0
#La fórmula para el modelo de crecimiento de Gompertz es la siguiente:
#N(t) = N0 * e^(-a * e^(-bt))
def GOMPERTZ(x,N0,a,b):
return N0 * np.exp(-a * np.exp(-b * x))
#La función logística se puede escribir de la siguiente manera:#
#N(t) = K / (1 + ae^(-rt))
def LOGISTICA(x,k,a,r):
return k / (1 + a * np.exp(-r * x))
def SIGMOID(x, L, k ,x0):
"""
Modelo sigmoidal (función logística).
Parámetros:
- x: Array de valores independientes.
- L: Límite superior de la curva.
- x0: Punto medio de la curva en el eje x.
- k: Pendiente de la curva en el punto medio.
Retorna:
- Array con los valores de la función sigmoidal para cada valor de x.
"""
return L / (1 + np.exp(-k * (x - x0)))
def GOMPERTZ(x,N0,a,b):
return N0 * np.exp(-a * np.exp(-b * x))
def GOMPERTZ2(x,N0,p,b):
return N0 * np.exp(- np.exp(p -b * x))
def relacionPolimorfica(x,a,b,c):
return a * ( 1 - np.exp( -b * x)) ** c
def MITSCHERLICH(x,a,L,b):
return a * ( 1 - L* np.exp( -b * x))
#
def ChapmanRichards(x,a,b,c):
return a * (1 -np.exp( -b * x)) ** c + np.exp(1)
def HossfeldI(x,a,b,c):
return ((x ** 2) / (a + b*x) ** 2) + np.exp(1)
def Schumacher(x,a,b,c):
return a * np.exp(-b * ( 1/ x)) + np.exp(1)
def Weibull(x,a,b,c):
return a * ( 1 - np.exp( -b * x ** c)) + np.exp(1)
En esta sección se definen un conjunto de funciones para encontrar el mejor modelo:
def evaluar_modelo(funcion, X, Y,popt):
predicciones = funcion(X, *popt)
# Calcular métricas básicas
mse = mean_squared_error(Y, predicciones)
r2 = r2_score(Y, predicciones)
mae = mean_absolute_error(Y, predicciones)
mape = np.mean(np.abs((Y - predicciones) / Y)) * 100
std_residuos = np.std(Y - predicciones)
# Otras métricas adicionales
suma_residuos_abs = np.sum(np.abs(Y - predicciones))
max_error_abs = np.max(np.abs(Y - predicciones))
mean_error_abs = np.mean(np.abs(Y - predicciones))
# Métricas adicionales
suma_residuos_cuad = np.sum((Y - predicciones)**2)
mediana_residuos = np.median(Y - predicciones)
percentil_90 = np.percentile(Y - predicciones, 90)
varianza_residuos = np.var(Y - predicciones)
suma_cuad_residual_abs = np.sum(np.abs((Y - predicciones)) ** 2)
min_error_abs = np.min(np.abs(Y - predicciones))
skew_residuos = scipy.stats.skew(Y - predicciones)
kurtosis_residuos = scipy.stats.kurtosis(Y - predicciones)
coef_variacion_residuos = np.std(Y - predicciones) / np.mean(Y - predicciones)
error_cuadratico_medio = np.sqrt(np.mean((Y - predicciones) ** 2))
metricas = {
"MSE": mse,
"R2": r2,
"MAE": mae,
"MAPE": mape,
"STD Residuos": std_residuos,
"Suma Residuos Abs": suma_residuos_abs,
"Max Error Abs": max_error_abs,
"Mean Error Abs": mean_error_abs,
"Suma Residuos Cuadráticos": suma_residuos_cuad,
"Mediana Residuos": mediana_residuos,
"Percentil 90": percentil_90,
"Varianza Residuos": varianza_residuos,
"Suma Cuadrada de los Residuos Absolutos": suma_cuad_residual_abs,
"Min Error Abs": min_error_abs,
"Skewness de los Residuos": skew_residuos,
"Kurtosis de los Residuos": kurtosis_residuos,
"Coeficiente de Variación de los Residuos": coef_variacion_residuos,
"Error Cuadrático Medio": error_cuadratico_medio
}
return metricas
def getMinimosCuadrados(funcion,X,Y):
def residuals(params, x, y):
return y - funcion(x, *params)
initial_guess = [1,1,1]
result = least_squares(residuals, initial_guess, args=(X, Y))
N0, t0, r = result.x
return result.x
# Establecer la visualización en línea y mejorar la calidad de la imagen
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
def getGrafico2(funcion, X, Y):
# Generar puntos para el ajuste
t_fit = np.linspace(min(X), max(X), 100)
# Obtener los parámetros del ajuste
a,b,c = getMinimosCuadrados(funcion,X,Y)
print(funcion.__name__,a,b,c)
# Evaluar la función ajustada en los puntos de ajuste
metricas = evaluar_modelo(funcion, X, Y,[a,b,c])
N_fit = funcion(t_fit,a,b,c)
# Plotear los datos y el ajuste
plt.figure(figsize=(15, 9)) # Tamaño de la figura
plt.scatter(X, Y, color='blue', label='Datos') # Plotear los datos
plt.plot(t_fit, N_fit, 'r-', label=f'Ajuste {funcion.__name__}') # Plotear el ajuste
plt.xlabel('Tiempo [Años]')
plt.ylabel('DAP [mm]')
plt.legend()
plt.grid(True) # Agregar una rejilla
plt.title('Gráfico de Ajuste') # Título del gráfico
plt.tight_layout() # Ajustar el diseño de la figura para evitar superposiciones
plt.show()
# Obtener métricas del modelo
#metricas = evaluar_modelo(funcion, X, Y, params)
metricas["Funcion"] = funcion.__name__
metricas["b0"] = a
metricas["b1"] = b
metricas["b2"] = c
return metricas
def MODELOS2(X,Y):
d1 = getGrafico2(SWEDA,X,Y)
d2 = getGrafico2(GOMPERTZ,X,Y)
d3 = getGrafico2(LOGISTICA,X,Y)
d4 = getGrafico2(SIGMOID,X,Y)
d5 = getGrafico2(GOMPERTZ2,X,Y)
d6 = getGrafico2(relacionPolimorfica,X,Y)
d7 = getGrafico2(MITSCHERLICH,X,Y)
d8 = getGrafico2(ChapmanRichards,X,Y)
d9 = getGrafico2(HossfeldI,X,Y)
d10 =getGrafico2(Schumacher,X,Y)
d11 =getGrafico2(Weibull,X,Y)
return [d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11]
dfClean = df[['Clasificacion Densidad','EDAD','DAP [mm]']].dropna()
Esta es la parte del proceso más larga y que requiere mayor capacidad de cómputo.
acumulador = []
for i in dfClean['Clasificacion Densidad'].sort_values().unique():
print(i)
aux2 = dfClean[dfClean['Clasificacion Densidad'] == i]
auxX = np.array(aux2['EDAD'])
auxY = np.array(aux2['DAP [mm]'])
data = pd.DataFrame(MODELOS2(auxX,auxY))
data["Clasificacion Densidad"] = i
acumulador.append(data.copy())
break
errores = pd.concat(acumulador)
errores["Fecha_Analisis"] = datetime.datetime.now()
1 SWEDA 21.04119438647293 0.046419999139448886 7.10587799699723
GOMPERTZ 38.46327179557505 1.725621174054583 0.13785611393744854
LOGISTICA 37.032921601250976 3.2109256850486454 0.18705747048150667
SIGMOID 37.03315033104063 0.18705474512879688 6.236430348041789
GOMPERTZ2 38.46292029061756 0.5455947007342655 0.13785926728915074
C:\Users\limc_\AppData\Local\Temp\ipykernel_25800\2818810145.py:37: RuntimeWarning: invalid value encountered in power return a * ( 1 - np.exp( -b * x)) ** c
relacionPolimorfica 44.160254594609384 0.0685756872552295 0.8134078716935539
MITSCHERLICH 41.366632283607444 0.9515908033103654 0.08787302045324784
C:\Users\limc_\AppData\Local\Temp\ipykernel_25800\2818810145.py:44: RuntimeWarning: invalid value encountered in power return a * (1 -np.exp( -b * x)) ** c + np.exp(1)
ChapmanRichards 40.27005479755286 0.07826878248230645 0.9718266320152252
HossfeldI -0.6570581261937979 -0.1451961919079237 1.0
Schumacher 41.20107176680019 5.904595978915142 1.0
Weibull 41.09688061447227 0.08410061768487219 0.9664273943593613
errores
| MSE | R2 | MAE | MAPE | STD Residuos | Suma Residuos Abs | Max Error Abs | Mean Error Abs | Suma Residuos Cuadráticos | Mediana Residuos | ... | Skewness de los Residuos | Kurtosis de los Residuos | Coeficiente de Variación de los Residuos | Error Cuadrático Medio | Funcion | b0 | b1 | b2 | Clasificacion Densidad | Fecha_Analisis | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.558134 | 0.468702 | 4.029327 | 17.917199 | 5.249462 | 12434.504221 | 25.503244 | 4.029327 | 85044.402026 | -0.300545 | ... | 0.368386 | 1.467820 | -1.467915e+02 | 5.249584 | SWEDA | 21.041194 | 0.046420 | 7.105878 | 1 | 2024-05-02 16:56:46.900983 |
| 1 | 24.920378 | 0.519556 | 3.803969 | 16.033188 | 4.992028 | 11739.048860 | 24.255991 | 3.803969 | 76904.285665 | -0.309014 | ... | 0.386523 | 1.732474 | -8.306736e+02 | 4.992031 | GOMPERTZ | 38.463272 | 1.725621 | 0.137856 | 1 | 2024-05-02 16:56:46.900983 |
| 2 | 25.099417 | 0.516104 | 3.824390 | 16.227436 | 5.009920 | 11802.066760 | 24.268368 | 3.824390 | 77456.799770 | -0.294159 | ... | 0.381045 | 1.697224 | -4.531408e+02 | 5.009932 | LOGISTICA | 37.032922 | 3.210926 | 0.187057 | 1 | 2024-05-02 16:56:46.900983 |
| 3 | 25.099417 | 0.516104 | 3.824391 | 16.227448 | 5.009920 | 11802.069612 | 24.268388 | 3.824391 | 77456.799764 | -0.294147 | ... | 0.381043 | 1.697223 | -4.531715e+02 | 5.009932 | SIGMOID | 37.033150 | 0.187055 | 6.236430 | 1 | 2024-05-02 16:56:46.900983 |
| 4 | 24.920378 | 0.519556 | 3.803968 | 16.033172 | 4.992028 | 11739.045338 | 24.255969 | 3.803968 | 76904.285673 | -0.308969 | ... | 0.386525 | 1.732475 | -8.306188e+02 | 4.992031 | GOMPERTZ2 | 38.462920 | 0.545595 | 0.137859 | 1 | 2024-05-02 16:56:46.900983 |
| 5 | 24.692280 | 0.523954 | 3.777398 | 15.808899 | 4.969129 | 11657.049170 | 24.326984 | 3.777398 | 76200.377269 | -0.262152 | ... | 0.398153 | 1.787827 | -8.413321e+02 | 4.969133 | relacionPolimorfica | 44.160255 | 0.068576 | 0.813408 | 1 | 2024-05-02 16:56:46.900983 |
| 6 | 24.727972 | 0.523266 | 3.779970 | 15.795186 | 4.972723 | 11664.988869 | 24.263190 | 3.779970 | 76310.521685 | -0.236278 | ... | 0.391794 | 1.774221 | 7.188269e+06 | 4.972723 | MITSCHERLICH | 41.366632 | 0.951591 | 0.087873 | 1 | 2024-05-02 16:56:46.900983 |
| 7 | 24.730752 | 0.523212 | 3.782651 | 15.864034 | 4.972995 | 11673.260247 | 24.324262 | 3.782651 | 76319.100526 | -0.287542 | ... | 0.397526 | 1.779102 | -5.878602e+02 | 4.973002 | ChapmanRichards | 40.270055 | 0.078269 | 0.971827 | 1 | 2024-05-02 16:56:46.900983 |
| 8 | 24.684556 | 0.524103 | 3.771243 | 15.669904 | 4.968355 | 11638.056309 | 24.195639 | 3.771243 | 76176.539028 | -0.198491 | ... | 0.401314 | 1.786289 | 2.991313e+03 | 4.968355 | HossfeldI | -0.657058 | -0.145196 | 1.000000 | 1 | 2024-05-02 16:56:46.900983 |
| 9 | 24.873829 | 0.520454 | 3.782467 | 15.595337 | 4.987348 | 11672.694631 | 24.035584 | 3.782467 | 76760.637141 | -0.173106 | ... | 0.413218 | 1.758948 | 3.638694e+02 | 4.987367 | Schumacher | 41.201072 | 5.904596 | 1.000000 | 1 | 2024-05-02 16:56:46.900983 |
| 10 | 24.728393 | 0.523257 | 3.782949 | 15.876494 | 4.972754 | 11674.181161 | 24.340027 | 3.782949 | 76311.819390 | -0.284092 | ... | 0.399331 | 1.781259 | -4.749267e+02 | 4.972765 | Weibull | 41.096881 | 0.084101 | 0.966427 | 1 | 2024-05-02 16:56:46.900983 |
11 rows × 24 columns
errores.to_sql("data_error", engine, if_exists='append', index=False)
55